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I. INTRODUCTION 



A basic building block for describing underwater explosions is the 
hydrodynamic properties across a shock wave in sea water. Underwater 
explosions were very extensively studied in World War II. Results from 
the World War II studies are reported in Underwater Explosion Research [1] . 
Richardson, Arons, and Halverson [2] discuss the calculation of hydro- 
dynamic properties of sea water at the front of a shock wave. Background 
information on the calculations appears in Chapter 2 of Cole [3]. 

Richardson, et al . , [2] used graphical techniques to interpolate 
thermodynamic data; the methods were crude and extremely tedious. With the 
advent of the modern programmable hand-held calculator, the results 
obtained by Richardson, et al . , [2] can be obtained with ease. This report 
discusses a program for the HP41CV, which uses the data of Richardson, 
et al . , [2] to calculate shock wave velocity, water velocity, specific 
volume, etc. Reference [2] is Appendix A. 

Since World War II, research has continued in the area of underwater 
explosions. A more recent survey by Professor Holt [4] appears in the 
Annual Reviews of Fluid Mechanics . Holt [4] discusses underwater nuclear 
explosions as well as chemical explosions. 

In passing, an interesting observation can be made concerning the level 
of talent working on underwater explosions in World War II. The names 
include G. I. Taylor, P. Bridgman, H. Bethe, J. von Neumann, and L. I. Sedov. 
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II. OUTLINE OF CALCULATION PROCEDURE 



A. Iteration of Properties with Rankine Hugoniot Relations 

The computer program has an iteration which is now discussed. The 
enthalpy jump across the shock wave is given by 

AH = -|(p - Pq) (v + Vq) (2.3) 

The same equation numbers and the same symbols used by Richardson, et al . , [2] 
are used here. Since p >> p^, p^ can be neglected. Also AH is given by 

AH =* a) + h (2.20) 

where oj is the undissipated enthalpy and h is the dissipated enthalpy. The 
dissipated enthalpy can be calculated two ways. First, 

a AH - oj 

Second, 

T„ 



A H =* h_ 
P T 



f T C dT 

L 0 



T, P 



(2.9) 



The symbol h^, which does not appear in reference [2], denotes the value 
calculated from AH - co. The symbol h^ is the value calculated from 
equation (2.9). Temperatures Tq and T^ are defined in reference [2]. The 
undissipated enthalpy is a function of specific volume, v, behind the shock 
wave; the function is 



0 ) 



2 n-1 

C 1 V 1 

- 1 ! 



(2.19) 



The Taic equation of state relating specific volume and pressure for water is 

p - Atsnoyv) 11 - i] 

where A(S] is a function only of entropy. Further 



A[S] » c^/nv^ 



(2.14) 



o 



Combining aquations (2.3), (2.19), and (2.14) yields 

2 2 
c c 

"h ■ :<V v,n - 11(v + V - - 11 <: - 21> 

Equation (2.21) indicates that h, is a function of v and v . The specific 

d 1 

volume, v^, is a function of T . By varying and v, the value of 
varies . 

When conditions across the shock wave are correct, 

w = v^v 

One varies v and until equality is achieved; since v = v(p,T^), only 
needs to be varied. 

This is a very brief description of the iteration. The reader should 
refer to reference [2] for more information. 

B. Temperature Behind Shock Wave 

The temperature increase behind the shock wave was calculated using 
equation (II-4) of reference [2]. Although the equation is somewhat 
lengthy, the calculation is straightforward except for two quantities: 
3'(tg) and 3 q . A cubic polynomial is given for B(tg); this polynomial was 
differentiated to obtain 



3'(t 0 ) « — ) 
d t 



t=t. 



The quantity 3 q is the coefficient of thermal expansion times specific 
volume and was evaluated using 




where v is specific volume. 



The symbol t is used for °C and T for °K. 
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III. DESCRIPTION OF COMPUTER PROGRAM 



A. Computer Program Flow Chart 

A flow chart of the program is given in Figure 1. An iteration loop 
3 HH 3 h^ 3 ET. The input "LOOP" is the criterion 
for acceptable difference between HH and HT. Since HH and HT have units 
of Joule/kg, "LOOP 11 also has units. The iteration variable is t, , which 
is the adiabatic temperature in °C. 

When an acceptable value for t^ has been found, the program auto- 

3 

matically goes to subroutine SHOCK. SHOCK calculates v, cm /gm; u, m/sec; 

U, m/sec; and c, m/sec. 

Following SHOCK, the program moves to subroutine DEL T, which calculates 
the difference in temperature from behind to front of shock wave. A variety 
of subroutines are used by DEL T including CP, 3TPRIME, Gl, and D2. 

3. Assignment of Storage Registers 

For the program, 34 storage registers are used. The assignment 
of variables to the registers is given in Table I. 

C. Program Listing 

Appendix 3 is a listing of the program. 

IV. OPERATION OF PROGRAM 

The program has been recorded on magnetic cards and can be inserted 
into the HP41CV using a HP82104A card reader. The inputs to the program 
will now be described. 

1. Turn on the HP41CV and activate USER. 

2. Assign the program to a key on the key board by pressing 

lUl ALPHA HH ALPHA, SIN 



occurs for making 



I INPUT LOOP ! 
I INPUT N=n ' 



CALCULATE 



1 


n+1 1 1 

n-1 ’ n’ a 


i 

i 






! 




INPUT INITIAL 
TEMPERATURE IN°( 


: 


i 










INPUT FINAL 
PRESSURE IN KILOBARS 






- 1 








CALCULATE v Q 








XEQ VW 





I 

INPUT t x 



XEQ Y 

CALCULATE (v /v) n 
FOR USE IN EQ. (2.21) 



XEQ VW 
CALCULATE v 



XEQ SPEED 
CALCULATE c^/n 



XEQ HT 
CALCULATE h T 



CALCULATE 



DISPLAY HH, HT 

AND HH-HT * 



I ABS (HH-HT) < LOOP~ 

r 

YES 



I ' GO TO "SHOCK" ! 

: CALCULATE SPECIFIC 
VOLUME v 3EHIND 
SHOCK WAVE 



j CALCULATE SPEED 
, OF WATER 3EHIND ] 

: SHOCK WAVE, u ' 

CALCULATE SPEED : 
OF SHOCK WAVE. U i 

| CALCULATE SPEED ' 

; OF SOUND 3EHIND j 
• SHOCK WAVE , c ! 



XEQ DEL T 

I 

1 

XEQ C? 

CALCULATE C 

2 _ 



XEQ BT PRIME 
CALCULATE dB/dt 



XEQ G1 
CALCULATE G 



XEQ D2 
CALCULATE D 



1 CALCULATE AT 

: 

| DISPLAY AT 1 
1 ~END~ | 



Figure 1. Computer Flow Chart 



Table I. Storage Registers 



Register Symbol 


Def inition 


Units 


Programs 


01 


C 1 


ad iaba tic t emper a tur e 


°C 


HT, VW, HH, 
SPEED 


02 


(c-55) 2 


quantity in polynomial 


°C 


3T 


03 


(t-25) 


quantity in polynomial 


°c 


VW 


04 


( t-55) 


quantity in polynomial 


°c 


3T 


05 


C 0 


temperature in front of 
shock wave 


°c 


HH 


06 


h T 


HT, dissipated enthalpy 


Joule/kg 


HT 


07 


bn 


HE, dissipated enthalpy 


Joule/kg 


HH 


08 


B(t) 


quantity in Tait equation 
of state 


kilobars 


BT 


09 


3 (t) 


quantity in Tait equation 
of state 


N/m 2 


3T 


10 


t-55 


quantity in polynomial 


°C 


BT PRIME 


11 


A£S] 


quantity in Tait equation 
of state 


N/m 2 


SPEED 


12 


P 


pressure behind shock wave 


kilobars 


Y 


13 


V 


specific volume 


cm^/gm 


SHOCK, VW 


14 


u 


velocity of water behind 
shock wave 


m/sec 


SHOCK 


15 


p 


pressure behind shock wave 


N/m^ 


SHOCK 


16 


U 


velocity of shock wave 


m/ sec 


SHOCK 


17 


C 1 


speed of sound at adiabatic 
temperature 


m/ sec 


SHOCK 


18 


c 


speed of sound behind shock 
wave 


m/ sec 


SHOCK 


19 


3 0 


coefficient of thermal 







expansion times specific 

volume v Q ta /kg°K BETA 
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Table 1 Continued. Storage Registers 



Register Symbol 


Definition 


Units 


Programs 


20 


2. 

c^/n 


quantitv in equation 
(2.21) 


2 2 
n“/sec" 


SPEED 


21 


(n+1) / (n-1) 


function of n 


- 


HH 


22 


1-1/n 


function of n 


- 


HH 


23 


-1/n 


function of n 


- 


HH 


24 


dB/dt 


derivative of B(t) 


N/m 2 °K 


3T PRIME 


25 


y 


(v^/v) n ratio of 
specific volumes 


- 


HH, Y 


26 


v i 


specific volume at 
ad ia ba t ic t emp er a t ur e 


cm^/gm 


HH, VW 


27 


v o 


specific volume in 
front of shock wave 


cm'Vgm 


HH 


28 


n 


exponent in Tait equa- 
tion of state 


- 


HH 


29 


C 

P 


heat capacity at 
constant pressure 


Joule/kg 

°K 


C? 


30 


G 


quantity in equation 
(II-4) 


°K 


G2, DEL T 


31 


D 


quantity in equation 
(II-4) 


- 


D2, DEL T 


32 


1+p/B 


quantity in Tait equa- 
tion of state 


- 


HH 


33 


(l+p/B) 1/n 


quantity in Tait equa- 
tion of state 


- 


HH 


34 


LOOP 


criterion for acceptable 
difference between HH 
and HT; typical value 
1.0 


Joule/kg 


HH 
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In this example, Che program has been assigned co the SIN key 
with blue H. 

3. Press SIN and observe LOOP?. LOOP is the value of HH-HT which 
is acceptable. When ABS (HH-HT) is less than LOOP, the program 
is out of the iteration loop. A typical value for LOOP 

could be 1.0 Joule/kg. Values for a sample case are given in 
parentheses. (1.0) 

4. Press R/S and observe N? . N is the exponent in Tait equation 

of state. Holt [4] uses 7.0. Richardson, et al ., [2] use 7.15. 
Insert a value for N. (7.15) 

5. Press R/S and observe TEMP 0?. This is the temperature in front 
of the shock wave in °C. Insert a value. (0) 

6. Press R/S and observe PRESSURE?. This is the pressure behind the 
shock wave in kilobars. Insert a value. (80) 

7. Press R/S and observe TEMP 1?. This is the adiabatic temperature 

discussed in reference [2]. Insert a value in °C. (180) 

8. Press R/S and observe crows foot moving across the display. In 
10 seconds, a value for HH appears. (520,245) 

9. Press R/S and observe value for HT. (724,255) 

10. Press R/S and observe value for HH-HT. (-204,010) 

Since the value for HH-HT exceeds LOOP, a new value for 
t^ =* TEMP 1 is needed. 

11. Press R/S and observe TEMP 1?. To decide what value of TEMP 1 
to insert, look at equation (2.9); h T increases as T 1 or TEMP 1 
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increases . 


Although not readily 


apparent 


from equation (2. 


decreases 


as TEMP 1 increases 


. For the 


example, a small 


cable can be 


made as follows: 






TEMP 1 


HH 


HT 


HH-HT 


130 


520,245 


724,255 


- 204,010 


160 


569,862 


642,846 


- 72,984 


149.2 


591,574 


598,988 


- 7,414 


148.0 


593,821 


594,120 


- 299 


147.75 


594,285 


593,106 


1,180 


147.949 


593,916 


593,913 


3 


147.9496 


593,915 


593,915 


- 0 



12. Press R/S and observe che following: 

v = 0.6576 (cm'Vgm) 

WATER U = 1,638 (m/sec) 

SHOCK U = 4,891 (m/sec) 

WAVE C = 6,264 (m/sec) 

13. Press R/S and observe the following: 

DEL T = 371 ,24°K 

Hence the temperature behind the shock wave is 
T » T 0 + iT =* 273.16 + 371.24 * 644. 4°K 
The calculation of shock properties has now been completed. To calculate 
values for another set of Cq and p, press SIN key with blue H. The 
computer is in USER mode. 
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V. SAMPLE RESULTS 



A comparison can be made between results calculated with the program 
and values tabulated by Richardson, et al . , [2]. Table II is a comparison. 
Most values are within 1/2 per cent. 

Since Richardson, et al . , [2] used graphical techniques for interpola- 
tion, the calculator values probably are more accurate. However, the 
thermodynamic data for water have been represented by polynomials. Errors 
in curve fitting are undoubtedly a few per cent. Hence, one would expect 
the values predicted here to be correct within a few per cent. 

VI. MAGNETIC CARDS 

Copies of the magnetic cards containing this program are available 

on request from the following address: 

Distinguished Professor Allen E. Fuhs 
Department of Aeronautics, Code 67 
Naval Postgraduate School 
Monterey, CA 93940 

(408)-646-2943 

AV-878-2948 

VII. CALCULATION OF WATER PROPERTIES 
DURING ISENTROPIC EXPANSION FROM SHOCKED CONDITIONS 

In order to calculate the complete flow field due to a shaped charge 
jet penetrating water, the calculation of temperature, internal energy, and 
enthalpy during is en tropic expansion is necessary. Appendix C discusses two 
subroutines which were developed after the main program was written. The 
subroutines calculate internal energy and enthalpy, subroutine E-3AR, and 
temperature, subroutine ISEN T. 
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Table II. Properties of Sea Water at a Shock Front* 
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***Vai U es In right-hand column are from Table IV of reference [2]. 
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APPENDIX A 



HYDRODYNAMIC PROPERTIES OF SEA WATER AT THE FRONT OF A SHOCK WAVE* 



*This is reference [2] reproduced from reference [1], 
is provided as a convenience to the reader. 



The technical paper 



Reprinted from The Journal of Chemical Physics, Vul. 15, No. 11. 735-794. November, 1947 

Printed in U. 3. A. 



Hydrodynamic Properties of Sea Water at the Front of a Shock Wave* 

J. M. Richardson** 

Baker Laboratory , Cornell University, Ithaca, .V. Y. 

AND 

A. B. Arons*** .and R. R. Halverson**** 

Underwater Explosives Research Laboratory, Woods Hole Oceanographic Institution, Woods Hole , Massachusetts 

(Received June 2, 1947) 

The Rankine-Hugoniot relations have been applied to appropriate equation-of-state data 
in order to calculate dae propagation velocity, particle velocity, enthalpy increment, Riemann 
function, etc. at shock fronts of various amplitudes in sea water. One set of tables provides 
values over a wide pressure range (up to about 80 kilobars) and is principally intended for use 
in conjunction with theories of propagation of shock waves originated by underwater ex- 
plosions. A second set of tables contains values which are closely spaced up to pressures of 14 
kilobars. These are calculated with somewhat greater precision and are intended for use in con- 
nection with experimental measurements of particle and propagation velocities, etc. 



L INTRODUCTION 

T has long been recognized that the velocity 
of propagation of sound waves of finite 
amplitude in a fluid medium is a function of the 
pressure in the wave. Lamb 1 ascribes the early 

* The work described in this report was performed 
under National Defense Research Committee Contracts 
OEMsr-121 with Cornell University and OEMsr-569 with 
the Woods Hole Oceanographic Institution. 

•• Present address: Bell Telephone Laboratories. Inc., 
Murray Hill, N. J. 

*** Present address: Department of Physics, Stevens 
Institute of Technology, Hoboken, N. J. 

**** Deceased. 

1 H. Lamb, Hydrodynamics (Cambridge University Press, 
London, 1932) 6th Ed., p. 481. 



development of the theory to independent in 
tigations of Eamshaw and Riemann. Qua 
tively this work indicated that, since the hi$ 
pressure portions of a wave travel with gre 
velocity, an arbitrarily-shaped pressure puls 
finite amplitude must, during propagation, a 
its shape in such a manner as to build up in 
shock front. By applying the laws of conse: 
tion of mass, energy, and momentum to 
transfer of matter across the shock fr 
Rankine and Hugoniot obtained a set of tl 
relations among the five variables: press 
density, particle velocity (u), shock fi 
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